## Run for a grid of carbon adders, separate for oil and for gas
adders.mesh = pracma::meshgrid(seq(0,90, by=10),
                               seq(0,60, by=5))
# adders.mesh = pracma::meshgrid(seq(0,30, by=30),
#                                seq(0,25, by=25))
adders.mesh = lapply(adders.mesh, t) # flip the matrices so that the first element (oil adders) 
# varies in the first dimension (rows)
# and the second element (gas adders) vary in the second dimension (cols)
names(adders.mesh) = c('oil_adder','gas_adder')
adder.dim = dim(adders.mesh[[1]])

sim.adder.grid.mesh = array(list(), dim=adder.dim, 
                            dimnames=list(paste0(as.character(adders.mesh$oil[,1]),'dollars'),
                                          paste0(as.character(adders.mesh$gas[1,]),'dollars')))
library(progress)
pb <- progress_bar$new(format = " Running meshgrid sims [:bar] :percent ETA= :eta",
                       total = prod(adder.dim)*2, clear = FALSE, width= 60)

st.loop = Sys.time()
for (e in 1:2) {
  if (e==1) load(working.data.loc)
  if (e==2) load(working.data.loc.high.e)
  
  for (i in 1:adder.dim[1]) {
    for (j in 1:adder.dim[2]) {
      # for (j in adder.dim[2]) {
      oil.adder = adders.mesh$oil_adder[i,j]
      gas.adder = adders.mesh$gas_adder[i,j]
      sim.adder.grid.mesh[[i,j]] = simulate.policy(carb.add.fed.sim.oil=oil.adder,
                                                   carb.add.fed.sim.gas=gas.adder) 
      # message(paste0(i, ' done of ', prod(adder.dim)*2,'. ',  Sys.time()))
      pb$tick()
    }
  }
  if (e==1) {
    sim.adder.grid.base.elasts = sim.adder.grid.mesh
    save('sim.adder.grid.base.elasts', file=working.data.dir%&%'/sim_adder_grid_mesh_base-elasts_'%&%today()%&%'.RData')
  }
  if (e==2) {
    sim.adder.grid.high.elasts = sim.adder.grid.mesh
    save('sim.adder.grid.high.elasts', file=working.data.dir%&%'/sim_adder_grid_mesh_high-elasts_'%&%today()%&%'.RData')
  }
}
ed.loop = Sys.time()
difftime(ed.loop, st.loop)
load(working.data.loc)
# load(working.data.dir%&%'/sim_adder_grid_mesh_base-elasts_2021-09-30.RData')
# load(working.data.dir%&%'/sim_adder_grid_mesh_high-elasts_2021-09-30.RData')

dimnames(sim.adder.grid.high.elasts)

# Check that these are still the opt
rev.opt.pair.b = gov.revenue(sim.adder.grid.base.elasts[['40dollars','10dollars']], agg=T)$tot.rev.summary
rev.opt.pair.h = gov.revenue(sim.adder.grid.high.elasts[['40dollars','10dollars']], agg=T)$tot.rev.summary

## Output results
deltas.adders.mesh.base = apply(sim.adder.grid.base.elasts,
                                1:2,
                                FUN=function(s) list(calc.deltas.for.mesh(sim.adder.grid.base.elasts[['0dollars','0dollars']], s[[1]])))
deltas.adders.mesh.high = apply(sim.adder.grid.high.elasts,
                                1:2,
                                FUN=function(s) list(calc.deltas.for.mesh(sim.adder.grid.high.elasts[['0dollars','0dollars']], s[[1]])))


mesh.outputs = c(colnames(sum.table)[1:7], 'Leakage-Oil','Leakage-Gas', colnames(sum.table)[8])
sum.table.mesh.base <- sum.table.mesh.high <- array(NA, dim=c(adder.dim[1],adder.dim[2], length(mesh.outputs)),
                                                    dimnames=c(dimnames(sim.adder.grid.base.elasts), 
                                                               list(mesh.outputs))) # 10 tables: one for each output
for (i in 1:adder.dim[1]) for (j in 1:adder.dim[2]) {
  sum.table.mesh.base[i,j,] = extract.vars.for.table.sep.leakage(deltas.adders.mesh.base[[i,j]][[1]],
                                                                 nm='',
                                                                 sim.base.temp = sim.adder.grid.base.elasts[[1,1]],
                                                                 sim.list.temp = sim.adder.grid.base.elasts[[i,j]])
  sum.table.mesh.high[i,j,] = extract.vars.for.table.sep.leakage(deltas.adders.mesh.high[[i,j]][[1]],
                                                                 nm='',
                                                                 sim.base.temp = sim.adder.grid.high.elasts[[1,1]],
                                                                 sim.list.temp = sim.adder.grid.high.elasts[[i,j]])
}
# write workbooks
output.sheet.names = mesh.outputs

library(openxlsx)
wb = createWorkbook('CarbonAdderMeshgrid')
# addWorksheet(wb, sheetName='Base')
for (o in 1:length(output.sheet.names)) {
  for (e in 1:2) { 
    if (e==1) x = sum.table.mesh.base[,,o]
    if (e==2) x = sum.table.mesh.high[,,o]
    rownames(x) = dimnames(sum.table.mesh.base)[[1]]
    name.temp = paste(output.sheet.names[o],ifelse(e==1,'-Base','-High'),'Elasts')
    addWorksheet(wb, sheetName=name.temp)
    writeData(wb, sheet=name.temp, x, rowNames = T)
  }
}
saveWorkbook(wb, output.dir%&%'/Carbon adder grid results/CarbonAdderMeshGrid_'%&%today()%&%'.xlsx', overwrite=TRUE)
rm(wb)

## Output table of revenue and emissions effects
dat <- data.table(sum.table.mesh.base[,,'Roy and CA Revenue'])
dat[, oil.adder := rownames(sum.table.mesh.base[,,'Roy and CA Revenue'])]
dat = melt(dat, id.var='oil.adder')
setnames(dat, c('oil.adder','gas.adder','rev'))

dat.2 <- data.table(sum.table.mesh.base[,,'Global emis'])
dat.2[, oil.adder := rownames(sum.table.mesh.base[,,'Roy and CA Revenue'])]
dat.2 = melt(dat.2, id.var='oil.adder')
setnames(dat.2, c('oil.adder','gas.adder','emis_base'))

dat.3 <- data.table(sum.table.mesh.high[,,'Global emis'])
dat.3[, oil.adder := rownames(sum.table.mesh.high[,,'Roy and CA Revenue'])]
dat.3 = melt(dat.3, id.var='oil.adder')
setnames(dat.3, c('oil.adder','gas.adder','emis_high'))

dat.4 <- data.table(sum.table.mesh.high[,,'Roy and CA Revenue'])
dat.4[, oil.adder := rownames(sum.table.mesh.high[,,'Roy and CA Revenue'])]
dat.4 = melt(dat.4, id.var='oil.adder')
setnames(dat.4, c('oil.adder','gas.adder','rev_high'))

dat = merge(dat, dat.2, by=c('oil.adder','gas.adder'))
dat = merge(dat, dat.3, by=c('oil.adder','gas.adder'))
dat = merge(dat, dat.4, by=c('oil.adder','gas.adder'))

dat[, oil.adder:=as.numeric(gsub('dollars','',oil.adder))]
dat[, gas.adder:=as.numeric(gsub('dollars','',gas.adder))]
dat = dat[order(oil.adder, gas.adder)]

dat[, oil.adder:='$'%&%oil.adder]
dat[, gas.adder:='$'%&%gas.adder]

dat = dat[as.numeric(gsub('\\$','',gas.adder))<=50]
dat.sign = sign(dat$rev)
ifelse(dat.sign<0,'-','')
dat.sign = sign(dat$rev_high)
ifelse(dat.sign<0,'-','')

# $40/ton for oil, $10/ton gas, is revenue-maximizing
dat[which.max(rev)]
dat[which.max(rev_high)]

lighten <- function(color, factor=1.4){
  col <- col2rgb(color)
  col <- col*factor
  col <- rgb(t(as.matrix(apply(col, 1, function(x) if (x > 255) 255 else x))), maxColorValue=255)
  col
}
palette()
## plot revenue data
ggplot(dat, aes(gas.adder, oil.adder)) +
  ggtitle("Change in Royalty & Carbon Adder Revenue, Base Elasticity Case (2020-2050 Average)") +
  geom_tile(aes(fill = rev)) + 
  scale_x_discrete(limits=unique(dat$gas.adder))+
  geom_text(aes(label = paste0(ifelse(rev<0,'-',''),'$',abs(round(rev, 1)))), family='serif') + 
  # scale_fill_gradient(low = rgb(235,211,103, maxColorValue = 255), high = 3, name=bquote(Delta~'Revenues ($b/yr)') ) +
  # scale_fill_gradient2(low = lighten('red', 1.2), mid='white', high = 3, name=bquote(Delta~'Revenues ($b/yr)') ) +
  scale_fill_gradient2(low = lighten(rgb(255,104,70, maxColorValue = 255), 0.9), mid='white', high = rgb(80,177,97,maxColorValue = 255), name=bquote(Delta~'Revenues ($b/yr)') ) +
  xlab('Gas Adder ($/ton)')+ylab('Oil Adder ($/ton)')+
  theme(text=element_text(family='serif'))
ggsave(filename = output.dir%&%'/Carbon adder grid results/revenue_heatmap_base-elast_'%&%today()%&%'.pdf',
       width=9.04, height=4.17, units='in')

ggplot(dat, aes(gas.adder, oil.adder)) +
  ggtitle("Change in Royalty & Carbon Adder Revenue, High Elasticity Case (2020-2050 Average)") +
  geom_tile(aes(fill = rev_high)) + 
  scale_x_discrete(limits=unique(dat$gas.adder))+
  geom_text(aes(label = paste0(ifelse(rev_high<0,'-',''),'$',abs(round(rev_high, 1)))), family='serif') + 
  # scale_fill_gradient(low = rgb(235,211,103, maxColorValue = 255), high = 3, name=bquote(Delta~'Revenues ($b/yr)') ) +
  # scale_fill_gradient2(low = lighten('red', 1.2), mid='white', high = 3, name=bquote(Delta~'Revenues ($b/yr)') ) +
  scale_fill_gradient2(low = lighten(rgb(255,104,70, maxColorValue = 255), 0.9), mid='white', high = rgb(80,177,97,maxColorValue = 255), name=bquote(Delta~'Revenues ($b/yr)') ) +
  xlab('Gas Adder ($/ton)')+ylab('Oil Adder ($/ton)')+
  theme(text=element_text(family='serif'))
ggsave(filename = output.dir%&%'/Carbon adder grid results/revenue_heatmap_high-elast_'%&%today()%&%'.pdf',
       width=9.04, height=4.17, units='in')

baseline.rev = gov.revenue(sim.base)$tot.rev.summary[year(date)>=2020, mean(total.revenue)]
dat[, gross.rev:= rev + baseline.rev]

ggplot(dat, aes(gas.adder, oil.adder)) +
  ggtitle("Total Royalty & Carbon Adder Revenue, Base Elasticity Case (2020-2050 Average)") +
  geom_tile(aes(fill = gross.rev)) + 
  scale_x_discrete(limits=unique(dat$gas.adder))+
  geom_text(aes(label = paste0(ifelse(gross.rev<0,'-',''),'$',abs(round(gross.rev, 1)))), family='serif') + 
  # scale_fill_gradient(low = rgb(235,211,103, maxColorValue = 255), high = 3, name=bquote(Delta~'Revenues ($b/yr)') ) +
  scale_fill_gradient(low = 'white', high = 3, name=bquote('Revenues ($b/yr)') ) +
  xlab('Gas Adder ($/ton)')+ylab('Oil Adder ($/ton)')+
  theme(text=element_text(family='serif'))
ggsave(filename = output.dir%&%'/Carbon adder grid results/gross_revenue_heatmap_'%&%today()%&%'.pdf',
       width=9.04, height=4.17, units='in')

## plot emissions data - base
ggplot(dat, aes(gas.adder, oil.adder)) +
  ggtitle("Change in Emissions, Base Elasticity Case (2020-2050 Average)") +
  geom_tile(aes(fill = emis_base)) + 
  scale_x_discrete(limits=unique(dat$gas.adder))+
  geom_text(aes(label = round(emis_base, 0)), family='serif') + 
  scale_fill_gradient(low = 3, high = rgb(235,211,103, maxColorValue = 255), name=bquote(Delta~'MMTCO'[2]~'e/yr') )+
  xlab('Gas Adder ($/ton)')+ylab('Oil Adder ($/ton)')+
  theme(text=element_text(family='serif'))
ggsave(filename = output.dir%&%'/Carbon adder grid results/emis_base-elast_heatmap_'%&%today()%&%'.pdf',
       width=9.04, height=4.17, units='in')
## plot emissions data - high
ggplot(dat, aes(gas.adder, oil.adder)) +
  ggtitle("Change in Emissions, High Elasticity Case (2020-2050 Average)") +
  geom_tile(aes(fill = emis_high)) + 
  scale_x_discrete(limits=unique(dat$gas.adder))+
  geom_text(aes(label = round(emis_high, 0)), family='serif') + 
  scale_fill_gradient(low = 3, high = rgb(235,211,103, maxColorValue = 255), name=bquote(Delta~'MMTCO'[2]~'e/yr') )+
  xlab('Gas Adder ($/ton)')+ylab('Oil Adder ($/ton)')+
  theme(text=element_text(family='serif'))
ggsave(filename = output.dir%&%'/Carbon adder grid results/emis_high-elast_heatmap_'%&%today()%&%'.pdf',
       width=9.04, height=4.17, units='in')

